Stan ODE notation beginning model

Observation Model

March 10 – March 30 parameters estimation

Postulated densities

knitr::opts_chunk$set(echo = FALSE)
library(ggplot2)
library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

library(rstan)

## Loading required package: StanHeaders

## rstan (Version 2.18.2, GitRev: 2e1f913d3ca3)

## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)

mu_ <- function(lambda, ylabel= "mu_h"){
    # We suppose mu_. as a exponential r.v.
    x <- seq(from = 0, to = 1, by=.001)
    mu_h_ <- dexp(x, rate = lambda)
    df <- data.frame(x=x, mu_h = mu_h_)
    mean_ <- 1/lambda
    mean_time <- 1/mean_ 
    text <- paste(mean_time,"days")
        p <- ggplot(data=df,
                aes(x = x,
                    y = mu_h)) + 
        xlab("x")+
        ylab(ylabel) +
        geom_line(colour = "#000000") +
        geom_area(alpha = 0.6, fill = "lightgray") + 
        geom_vline(xintercept = mean_, linetype="dotted", 
                color = "red", size=0.5)+
        annotate(geom = "text", x = 1.30 * mean_, y=0, label = text,
              color="black")
    return(p)
}
#
gamma_ <- function(alpha, beta, ylabel="gamma_S"){
    # We suppose mu_. as a exponential r.v.
    theta <- 1 / beta 
    x <- seq(from = 0, to = 1, by=.001)
    gamma_s_ <- dgamma(x, shape = alpha, 
                        scale= theta)
    df <- data.frame(x=x, gamma_s = gamma_s_)
    mean_ <- alpha * theta
    mean_time <- 1/mean_ 
    text <- paste(mean_time,"days")
    p <- ggplot(data=df,
                aes(x = x,
                    y = gamma_s)) + 
        geom_line(colour = "#000000") +
        geom_area(alpha = 0.6, fill = "lightgray") +
        geom_vline(xintercept = mean_, linetype="dotted", 
                color = "red", size=0.5) + 
        annotate(geom="text", x=1.30 * mean_, y=0, label = text,
              color="black")
        return(p)
}
pi_mu_h <- mu_(lambda = 25, ylabel= "mu_h")
pi_mu_i_s <- mu_(lambda = 35, ylabel= "mu_i_s")
pi_gamma_s <- gamma_(alpha = 10, beta = 100,  ylabel= "gamma_S")
pi_gamma_a <- gamma_(alpha = 10, beta = 50,  ylabel= "gamma_A")
pi_gamma_h <- gamma_(alpha = 10, beta = 250,  ylabel= "gamma_H")

fig6 <- ggplotly(pi_mu_h)
fig7 <- ggplotly(pi_mu_i_s)
fig8 <- ggplotly(pi_gamma_s)
fig9 <- ggplotly(pi_gamma_a)
fig10 <- ggplotly(pi_gamma_h)
fig6

fig7

fig8

fig9

fig10

Stan estimation

Lockdown-Vaccination-Treatment-Model

Model from March 24 July 14

March 24 –Jun 01 parameters estimation

Stan ODE notation

Stan MODEL

## functions {
##     real[] seirvt(real t,  // time
##              real[] y,
##              real[] theta,
##              real[] x_r,
##              int[] x_i) {
## //
##     real dy_dt[8];
##     real s = y[1];
##     real e = y[2];
##     real i_s = y[3];
##     real i_a = y[4];
##     real h = y[5];
##     real r = y[6];
##     real d = y[7];
##     real n_star;
## //
##     real beta_s = theta[1];
##     real beta_a = theta[2];
##     real kappa = theta[3];
##     real p = theta[4];
##     real delta_h = theta[5];
##     real mu_i_s = theta[6];
##     real mu_h = theta[7];
##     real gamma_s = theta[8];
##     real gamma_a = theta[9];
##     real gamma_h = theta[10];
## //
##     real force_infection;
##     real mu = 3.913894e-05;
##     n_star = s + e + i_s + i_a + h + r;
## 
##     force_infection = (beta_s * i_s + beta_a * i_a) / n_star;
##     dy_dt[1] = mu * n_star - force_infection * s - mu * s;
##     dy_dt[2] = force_infection * s - (mu + kappa) * e;
##     dy_dt[3] = p * kappa * e - (gamma_s + mu_i_s + delta_h) * i_s;
##     dy_dt[4] = (1.0 - p) * kappa * e - (gamma_a + mu) * i_a;
##     dy_dt[5] = delta_h * i_s  - (gamma_h + mu_h + mu) * h;
##     dy_dt[6] = gamma_s * i_s + gamma_a * i_a  + gamma_h * h - mu * r;
##     dy_dt[7] = mu_i_s * i_s + mu_h * h;
##     dy_dt[8] = p * kappa * i_s;
##     return dy_dt;
##   }
## }
## data {
##     int<lower = 1> n_obs;       // number of days observed
##     int<lower = 1> n_theta;     // number of model parameters
##     int<lower = 1> n_difeq;     // number of differential equations
##     int<lower = 1> n_pop;       // population
##     int y[n_obs];               // data, total number of infected
##     real t0;                    // initial time point (zero)
##     real ts[n_obs];             // time points observed
## }
## 
## transformed data {
##   real x_r[0];
##   int x_i[0];
## }
## parameters {
##     real <lower = 5e-6> theta[n_theta];
##     // real <lower = 0, upper = 1> E0;
##     // real <lower = 0, upper=1> IA0;
##     real <lower = 0, upper = n_pop> E0;
##     real <lower = 0, upper = n_pop> IA0;
## }
## transformed parameters{
##     real y_hat[n_obs, n_difeq]; // solution from the ODE solver
##     real y_init[n_difeq];
##     // initial conditions
##     // real IS0 = 8.297217e-06;
##     real IS0 = 22.0;
##     real H0 = 0.0;
##     real R0 = 0.0;
##     real D0 = 0.0;
##     real CIS0 = 74.0;
##     y_init[1] = n_pop - (IS0 + IA0 + E0);
##     y_init[2] = E0;
##     y_init[3] = IS0;
##     y_init[4] = IA0;
##     y_init[5] = H0;
##     y_init[6] = R0;
##     y_init[7] = D0;
##     y_init[8] = CIS0;
##     y_hat = integrate_ode_rk45(seirvt, y_init, t0, ts, theta, x_r, x_i);
## }
## model {
##     real lambda[n_obs];             // poisson parameter
##                                     // priors
##     theta[1] ~ normal(1.0, 0.1);    // beta_s
##     theta[2] ~ normal(1.0, 0.1);    // beta_a
##     theta[3] ~ gamma(10, 40);       // kappa
##     theta[4] ~ uniform(0.1, 0.5);   // p
##     theta[5] ~ gamma(10, 40);       // delta_h
##     theta[6] ~ exponential(33.33333);   // mu_i_s
##     theta[7] ~ exponential(25);   // mu_h
##     theta[8] ~ gamma(10, 100);      // gamma_s
##     theta[9] ~ gamma(10, 50);       // gamma_a
##     theta[10] ~ gamma(10, 250);      // gamma_h
## //
##     // E0 ~ uniform(1.121246e-07, 3.363737e-05);
##     // IA0 ~ uniform(7.848719e-06, 3.363737e-05);
##     E0 ~ uniform(74, 2100);
##     IA0 ~ uniform(74, 2100);
## //    likelihood
##     for (i in 1:n_obs){
##        //lambda[i] = y_hat[i, 8] * n_pop;
##         lambda[i] = y_hat[i, 8];
##     }
##     y ~ poisson(lambda);
## }
## generated quantities {
##     real R_0;      // Basic reproduction number
##     real mu = 3.913894e-05;
##     real beta_s = theta[1];
##     real beta_a = theta[2];
##     real p = theta[4];
##     real gamma_s = theta[4];
##     real gamma_a = theta[5];
##     real kappa = theta[6];
##     R_0 = p * beta_s * kappa / ((gamma_s + mu) * (kappa + mu))
##         + (1.0 - p) * beta_a * kappa / ((gamma_a + mu) * (kappa + mu));
## }

Observation Model